
#drafted by Lauren Hodgson ( lhodgson86@gmail.com ...  )
#GNU General Public License .. feel free to use / distribute ... no warranties

#use for changing the extent of a small number of ascii files
#note: contains test image script
################################################################################

library(SDMTools)

wd='/home/ctbccr/Ramona/Substropo_data/';setwd(wd) #define your working directory; set working directory

base.asc = read.asc.gz('/home/jc165798/Climate/CIAS/Australia/1km/baseline.76to05/base.asc.gz') #read in the base asc file
pos = read.csv('/home/jc165798/Climate/CIAS/Australia/1km/baseline.76to05/base.positions.csv',as.is=TRUE) #read in the position files

#read in layers
files=list.files()
i=5
columns=NULL
for (tfile in files) {cat(tfile,'\n')
	
	tasc=read.asc(tfile)
	tname=strsplit(tfile,'\\.')[[1]][1]

	pos[,i] = extract.data(cbind(pos$lon,pos$lat),tasc)   #extract the data from the ascii and append to pos file

	#create new asciis
	newasc=base.asc; newasc[cbind(pos$row,pos$col)]=pos[,i] #bind pos data to base.asc
	write.asc(newasc,paste(tname,'.new.asc',sep=''))
	i=i+1
	columns=c(columns,tname)
}
colnames(pos[,5:ncol(pos)])=columns


###################
#test image script
old=read.asc('bdensity.asc') #read in old ascii
new=read.asc('bdensity.new.asc') #read in new ascii

cols = colorRampPalette(c('skyblue','slateblue','forestgreen','yellow','red'))(101)
pnts=cbind(x=c(112,116,116,112), y=c(-11,-11,-18.5,-18.5))
png('new.png', width=7, height=7, units='cm', res=300, pointsize=5, bg='white')
	image(new,ann=FALSE,axes=FALSE,col=cols)
	legend.gradient(pnts,cols=cols,limits=round(range(new,na.rm=T)), title='limits', cex=1)
	dev.off()
png('old.png', width=7, height=7, units='cm', res=300, pointsize=5, bg='white')
	image(old,ann=FALSE,axes=FALSE,col=cols)
	legend.gradient(pnts,cols=cols,limits=round(range(old,na.rm=T)), title='limits', cex=1)
	dev.off()
png('base.png', width=7, height=7, units='cm', res=300, pointsize=5, bg='white')
	image(base.asc,ann=FALSE,axes=FALSE,col=cols)
	legend.gradient(pnts,cols=cols,limits=round(range(base.asc,na.rm=T)), title='limits', cex=1)
	dev.off()
##################

	
	
	